#include "yocto/cliff/wksp1d.hpp"
#include "yocto/lua/lua-state.hpp"
#include "yocto/lua/lua-config.hpp"
#include "yocto/string/vfs-utils.hpp"
#include "yocto/math/types.hpp"
#include "yocto/math/ode/drvck.hpp"
#include "yocto/math/ode/stiff-drvrs.hpp"
#include "yocto/ios/ocstream.hpp"
#include "yocto/auto-ptr.hpp"
#include "yocto/wtime.hpp"
#include "yocto/string/conv.hpp"
#include "yocto/eta.hpp"
#include "yocto/z++/gzstream.hpp"

using namespace yocto;
using namespace cliff;
using namespace math;


extern void compare_db();

////////////////////////////////////////////////////////////////////////////////
//
// Types
//
////////////////////////////////////////////////////////////////////////////////
typedef double                 Real;
typedef coord1D                Coord;
typedef layout1D               Layout;
typedef region1D<Real>::type   Region;
typedef array1D<Real>          Array;
typedef vertex1D<Real>::type   Vertex;
typedef wksp1D<Real,Real>      Workspace; 
typedef ghosts_setup<Coord>    Ghosts;
typedef ode::drvck<Real>::type ODEint;
typedef ode::stiff_drvrs<Real>::type StiffODEint;
typedef ode::field<Real>::type DiffEq;
typedef ode::field<Real>::diff Jacobian;

static Real   Kw       = 0;
static Real   Dh       = 0;
static Real   Dw       = 0;
static Real   pH_left  = 0;
static Real   pH_right = 0;
static Real   k2       = 0;
static Real   kd       = 0;
static Real   dt0      = 1e-6;
static Real   tol      = 1e-3;
static wtime *chrono = NULL;

#define GFX 1
#if GFX == 1
#include <FL/gl.h>
#include <FL/FL_Gl_Window.h>
#include <FL/Fl.h>

namespace {
    
    class MyWindow : public  Fl_Gl_Window
    {
    public:
        const Workspace &wksp;
        MyWindow(int X, int Y, int W, int H, const char *L, const Workspace &wk)
        : Fl_Gl_Window(X, Y, W, H, L),
        wksp(wk)
        {
            std::cerr << "** Creating OpenGL Window" << std::endl;
        }
        
        virtual ~MyWindow() throw()
        {
            std::cerr << "** Destroy OpenGL Window" << std::endl;
        }
        
        
        virtual void draw()
        {
            if (!valid()) 
            {
                valid(1);
                glLoadIdentity();
                glViewport(0,0,w(),h());
            }
            
            glClear(GL_COLOR_BUFFER_BIT);
            display( wksp.X, wksp["w"], FL_GREEN  );
            display( wksp.X, wksp["h"], FL_RED    );
        }
        
        inline void display( const Array &X, const Array &Y, Fl_Color c )
        {
            Real xmin = wksp.region.vmin;
            Real xmax = wksp.region.vmax;
            Real vmin=0,vmax=0;
            Y.get_min_max(vmin, NULL, vmax, NULL);
            Real ymin = -Log10(vmax);
            Real ymax = -Log10(vmin);
            const Real xlen = xmax-xmin;
            const Real dx   = 0.05 * xlen;
            const Real ylen = ymax - ymin;
            const Real dy   = 0.05 * ylen;
            glPushMatrix();
            glLoadIdentity();
            glOrtho(xmin-dx,xmax+dx, ymin-dy, ymax+dy, 1, -1);
            gl_color(c);
            glBegin(GL_LINE_STRIP);
            for( unit_t i=wksp.lower; i <= wksp.upper; ++i )
            {
                glVertex3d(X[i], -Log10(Y[i]), 0);
            }
            glEnd();
            
            glPopMatrix();
        }
        
        
    private:
        YOCTO_DISABLE_COPY_AND_ASSIGN(MyWindow);
    };
    
}
#endif



////////////////////////////////////////////////////////////////////////////////
//
// initialize
//
////////////////////////////////////////////////////////////////////////////////
static inline void initialize( Workspace &W )
{
    const size_t volumes = W.upper;
    Array &h    = W["h"];
    Array &w    = W["w"];
    Array &step = W["step"];
    {
        h[0]    = Pow(10.0,-pH_left);
        w[0]    = Kw/h[0];
        step[0] = dt0;
        const Real h_right = Pow(10.0,-pH_right);
        const Real w_right = Kw/h_right;
        for( unit_t i=1; i <= volumes; ++i )
        {
            h[i]    = h_right;
            w[i]    = w_right;
            step[i] = dt0;
        }
    }
    
}


////////////////////////////////////////////////////////////////////////////////
//
// standard scheme
//
////////////////////////////////////////////////////////////////////////////////

//==============================================================================
//
// Differential Fields computation
//
//==============================================================================
static inline void compute_fields(const Array &X, 
                                  Array       &Flux,
                                  Array       &delta, 
                                  const Real   D,
                                  const Real   idx,
                                  const size_t volumes)
{
    const Real   D_over_dx = D * idx;
    const Real   fac = D_over_dx * 0.5;
    
    //-- bulk
    for( size_t i=1; i < volumes; ++i )
        Flux[i] = fac * ( X[i+1] - X[i-1] );
    
    //-- boundaries
    Flux[0]      =  fac * ( 4.0 * X[1]         - (X[2]        +3.0*X[0]      ) );
    Flux[volumes]= -fac * ( 4.0 * X[volumes-1] - (X[volumes-2]+3.0*X[volumes]) );
    
    //-- increases
    for( size_t i=0; i < volumes; ++i )
    {
        delta[i] = idx * (Flux[i+1] - Flux[i]);
    }
}

//==============================================================================
//
// explicit standard diffusion
//
//==============================================================================
static inline double explicit_standard_diffusion( Workspace &W, const Real dt)
{
    assert(chrono);
    const double t1 = chrono->query();
    
    Array       &h       =  W["h"];
    Array       &w       =  W["w"];
    Array       &Fh      =  W["Fh"];
    Array       &Fw      =  W["Fw"];
    Array       &delta_h =  W["delta_h"];
    Array       &delta_w =  W["delta_w"];
    const size_t volumes =  W.upper; assert(0==W.lower);
    const Real   half_dt = 0.5 * dt;
    
    //==========================================================================
    // compute differential fields
    //==========================================================================
    compute_fields( h, Fh, delta_h, Dh, W.inv_d, volumes );
    compute_fields( w, Fw, delta_w, Dw, W.inv_d, volumes );
    
    
    //==========================================================================
    // update concentrations
    //==========================================================================
    for( size_t i=1; i < volumes; ++i )
    {
        h[i] += half_dt * (delta_h[i-1]+delta_h[i]);
        w[i] += half_dt * (delta_w[i-1]+delta_w[i]);
    }
    
    return chrono->query() - t1;
}

//==============================================================================
//
// reactive equation
//
//==============================================================================
struct reactions
{
    void diffeq( array<Real> &dYdt, Real t, const array<Real> &Y )
    {
        const Real rate = (kd - k2 * Y[1] * Y[2]);
        dYdt[1] = rate; //! F[1]
        dYdt[2] = rate; //! F[2]
    }
    
    void jacobn( array<Real> &dFdt, matrix<Real> &dFdY, Real t, const array<Real> &Y )
    {
        dFdt[1] = 0;
        dFdt[2] = 0;
        
        dFdY[1][1] = dFdY[2][1] = -k2 * Y[2];
        dFdY[1][2] = dFdY[2][2] = -k2 * Y[1];
        
    }
};




//==============================================================================
//
// reactive part
//
//==============================================================================
static inline double reaction(Workspace   &W,
                              const Real   t, 
                              const Real   dt,
                              ODEint      &odeint,
                              DiffEq      &drvs,
                              array<Real> &Y )
{
    assert(chrono);
    const double t1 = chrono->query();
    const size_t volumes = W.upper; assert(0==W.lower);
    Array       &h       = W["h"];
    Array       &w       = W["w"];
    Array       &step    = W["step"];
    for( size_t i=1; i < volumes; ++i )
    {
        Real &h_i = h[i];
        Real &w_i = w[i];
        Y[1] = h_i;
        Y[2] = w_i;
        odeint(drvs,Y,t,t+dt,step[i]);
        h_i = Y[1];
        w_i = Y[2];
    }
    return chrono->query() - t1;
}

static inline double stiff_reaction(Workspace   &W,
                                    const Real   t, 
                                    const Real   dt,
                                    StiffODEint  &stiff_odeint,
                                    DiffEq       &derivs,
                                    Jacobian     &jacobn,
                                    array<Real>  &Y )
{
    assert(chrono);
    const double t1 = chrono->query();
    const size_t volumes = W.upper; assert(0==W.lower);
    Array       &h       = W["h"];
    Array       &w       = W["w"];
    Array       &step    = W["step"];
    for( size_t i=1; i < volumes; ++i )
    {
        Real &h_i = h[i];
        Real &w_i = w[i];
        Y[1] = h_i;
        Y[2] = w_i;
        stiff_odeint(derivs,jacobn,Y,t,t+dt,step[i]);
        h_i = Y[1];
        w_i = Y[2];
    }
    return chrono->query() - t1;
}


static inline Real compute_error( const Real h, const Real w )
{
    return Fabs( (h*w-Kw)/Kw );
}

static inline Real find_max_error( const Workspace &W )
{
    const size_t volumes = W.upper; assert(0==W.lower);
    const Array       &h       = W["h"];
    const Array       &w       = W["w"];
    Real ans = compute_error( h[0], w[0] );
    for( size_t i=1; i <= volumes; ++i )
    {
        const Real tmp = compute_error( h[i], w[i] );
        if( tmp > ans ) ans = tmp;
    }
    return ans;
}

////////////////////////////////////////////////////////////////////////////////
//
// constrained scheme
//
////////////////////////////////////////////////////////////////////////////////

//==============================================================================
// Newton Step
//==============================================================================
static size_t num_newton = 0;
static inline void newton( Real &h, Real &w )
{
    assert(h>0);
    assert(w>0);
    while( true )
    {
        ++num_newton;
        const Real Gam  = Kw - (h*w);
        const Real dxi  = Gam/(w+h);
        const Real dh   = dxi;
        const Real dw   = dxi;
        h += dh;
        w += dw;
        if( Fabs(dh) <= Fabs( tol * h ) && Fabs(dw) <= Fabs( tol * w ) )
            break;
    }
}

static inline void test_newton()
{
    Real h = 1;
    Real w = 0.2;
    newton( h, w );
    std::cerr << "h=" << h << ", w=" << w << ", K=" << w*h << std::endl;
}

static inline void try_explicit_constrained_diffusion( Workspace &W, const Real dt)
{
    
    Array       &h       =  W["h"];
    Array       &w       =  W["w"];
    Array       &Fh      =  W["Fh"];
    Array       &Fw      =  W["Fw"];
    Array       &delta_h =  W["delta_h"];
    Array       &delta_w =  W["delta_w"];
    Array       &cdh     =  W["cdh"];
    Array       &cdw     =  W["cdw"];                    
    
    const size_t volumes =  W.upper; assert(0==W.lower);
    
    const Real   half_dt = 0.5 * dt;
    
    //==========================================================================
    // compute differential fields
    //==========================================================================
    compute_fields( h, Fh, delta_h, Dh, W.inv_d, volumes );
    compute_fields( w, Fw, delta_w, Dw, W.inv_d, volumes );
    
    
    
    //==========================================================================
    // update concentrations
    //==========================================================================
    for( size_t i=1; i < volumes; ++i )
    {
        const Real dh  = half_dt * (delta_h[i-1]+delta_h[i]);
        const Real dw  = half_dt * (delta_w[i-1]+delta_w[i]);
        const Real h_i = h[i];
        const Real w_i = w[i];
        const Real fac = 1.0/(h_i+w_i);
        
        //-- constrained global increase
        cdh[i] = fac * h_i * (dh - dw);
        cdw[i] = fac * w_i * (dw - dh); 
        
        
    }
    
    
}

static inline bool is_valid_constrained_diffusion( const Workspace &W )
{
    const Array       &h       =  W["h"];
    const Array       &w       =  W["w"];
    const Array       &cdh     =  W["cdh"];
    const Array       &cdw     =  W["cdw"];    
    const size_t       volumes =  W.upper; assert(0==W.lower);
    
    for( size_t i=1; i < volumes; ++i )
    {
        if( cdh[i] <= -h[i] || cdw[i] <= -w[i] )
            return false;
    }
    return true;
    
}

static inline double explicit_constrained_diffusion( Workspace &W, const Real dt)
{
    assert(chrono);
    const double t1 = chrono->query();
    Array       &h       =  W["h"];
    Array       &w       =  W["w"];
    Array       &cdh     =  W["cdh"];
    Array       &cdw     =  W["cdw"];   
    const size_t volumes =  W.upper; assert(0==W.lower);
    
    Real dt_left = dt;
    while( dt_left > 0 )
    {
        Real dt_done = dt_left;
        try_explicit_constrained_diffusion(W,dt_done);
        while( !is_valid_constrained_diffusion(W) )
        {
            dt_done *= 0.5;
            for( size_t i=1; i < volumes; ++i )
            {
                cdh[i] *= 0.5;
                cdw[i] *= 0.5;
            }
        }
        
        for( size_t i=1; i < volumes; ++i )
        {
            h[i] += cdh[i]; assert( h[i] > 0 );
            w[i] += cdw[i]; assert( w[i] > 0 );
            newton( h[i], w[i] );
        }
        
        dt_left -= dt_done;
    }
    return chrono->query() - t1;
}


////////////////////////////////////////////////////////////////////////////////
//
// output
//
////////////////////////////////////////////////////////////////////////////////
static inline void save_fields( const string &filename, const Workspace &W )
{
    ios::ocstream fp(filename,false);
    for( unit_t i=W.lower;i<=W.upper;++i)
    {
        fp("%g %g %g %g %g\n",W.X[i], W["h"][i],W["w"][i],W["Fh"][i],W["Fw"][i]);
    }
}


static inline void show_time( const double t, const eta &progress, const double rate )
{
    fprintf( stderr, "t= %10.5f  : ETA: %10.2f s   | fps= %7.1f    \r", t, progress.time_left, rate );fflush(stderr);
}



static inline void save_db(const string    & filename, 
                           const double      stamp,
                           const Workspace & W, 
                           bool              append,
                           const double      t_cycle)
{
    ios::ogzstream fp( filename, 9, append );
    fp.emit<double>( stamp );
    fp.emit<uint32_t>( W.width );
    fp.save( &W["h"][0], W.width * sizeof(Real) );
    fp.emit<double>( t_cycle );
}

////////////////////////////////////////////////////////////////////////////////
//
// main code
//
////////////////////////////////////////////////////////////////////////////////

int main( int argc, char *argv[] )
{
    const char *prog = _vfs::get_base_name(argv[0]);
    try 
    {
        
        wtime __chrono;
        chrono = & __chrono;
        chrono->start();
        
        ////////////////////////////////////////////////////////////////////////
        //
        // Load config file
        //
        ////////////////////////////////////////////////////////////////////////
        std::cerr << "-- Reading Configuration" << std::endl;
        if( argc <= 4) 
            throw exception("usage: %s config.lua dt_save tmax tol", prog );
        Lua::State VM;
        lua_State *L = VM();
        Lua::Config::DoFile(L,argv[1]);
        Real dt_save = strconv::to_real<Real>( argv[2], "dt_save"); 
        Real tmax    = strconv::to_real<Real>( argv[3], "tmax" );
        tol          = strconv::to_real<Real>( argv[4], "tol" );
        
        ////////////////////////////////////////////////////////////////////////
        //
        // Fetch Parameters
        //
        ////////////////////////////////////////////////////////////////////////
        std::cerr << "-- Reading Parameters" << std::endl;
        const Region sim_region(0.0,Real(  Lua::Config::Get<lua_Number>(L,"length")  ));
        const Ghosts sim_ghosts;
        const Layout sim_layout(0,  Coord( Lua::Config::Get<lua_Number>(L,"volumes") ));
        Kw       = Real(  Lua::Config::Get<lua_Number>(L,"Kw")       );
        Dh       = Real(  Lua::Config::Get<lua_Number>(L,"Dh")       );
        Dw       = Real(  Lua::Config::Get<lua_Number>(L,"Dw")       );
        pH_left  = Real(  Lua::Config::Get<lua_Number>(L,"pH_left")  );
        pH_right = Real(  Lua::Config::Get<lua_Number>(L,"pH_right") );
        kd       = Real(  Lua::Config::Get<lua_Number>(L,"kd")       );
        k2       = Real(  Lua::Config::Get<lua_Number>(L,"k2")       );
        const string react = Lua::Config::Get<string>( L, "reaction" );
        bool  explicit_reaction = true;
        if( react != "explicit" )
        {
            if( react == "implicit" )
            {
                explicit_reaction = false;
            }
            else {
                throw exception("Invalid reaction='%s'", &react[0] );
            }
        }
        std::cerr << "-- reaction: " << (explicit_reaction ? "explicit" : "implicit") << std::endl;
        const size_t volumes = sim_layout.upper; assert(0==sim_layout.lower);
        
        ////////////////////////////////////////////////////////////////////////
        //
        // Create workspace: for species, fluxe, increases,
        // plus one for the control step
        //
        ////////////////////////////////////////////////////////////////////////
        std::cerr << "-- Creating Workspace with " << volumes << " volumes" << std::endl;
        const char   *var_names[] = { "h", "w", "Fh", "Fw", "delta_h", "delta_w", "cdh", "cdw", "step" };
        const size_t  var_count   = sizeof(var_names)/sizeof(var_names[0]);
        Workspace     W( sim_layout, sim_ghosts, sim_region , var_count, var_names );
        
        if( false )
        {
            std::cerr << "-- Testing Newton..." << std::endl;
            test_newton();
            return 0;
        }
        
        ////////////////////////////////////////////////////////////////////////
        //
        // Compute dt
        //
        ////////////////////////////////////////////////////////////////////////
        const Real length_sq      =  sim_region.length * sim_region.length;
        const Real fast_diff_time =  length_sq / Dh;
        const Real slow_diff_time =  length_sq / Dw;
        std::cerr << "-- FastDiffTime= " << fast_diff_time << std::endl;
        std::cerr << "-- SlowDiffTime= " << slow_diff_time << std::endl;
        const Real dt_max = 0.1 * W.delta * W.delta / Dh;
        std::cerr << "-- dt_max      = " << dt_max << std::endl;
        const Real dt_rnd = Pow(10.0,Floor(Log10(dt_max)));
        const Real dt_fac = Floor( dt_max/dt_rnd );
        const Real dt     = dt_fac * dt_rnd;
        std::cerr << "-- dt          = " << dt << std::endl;
        
        
#if GFX == 1
        const int MySize = 512;
        auto_ptr<MyWindow> win( new MyWindow(10,10,MySize,MySize,"Simulation",W) );
        win->mode( FL_RGB );
        win->show();
        Fl::check();
#endif
        
        ////////////////////////////////////////////////////////////////////////
        //
        // Adaptive solver and differential system
        //
        ////////////////////////////////////////////////////////////////////////
        const size_t nvar   = 2;
        ODEint       odeint;
        StiffODEint  stiff_odeint;
        odeint.eps          = tol;
        stiff_odeint.eps    = tol;
        odeint.start(nvar);
        stiff_odeint.start(nvar);
        
        vector<Real> Y(nvar,0.0);
        reactions    rxn;
        DiffEq       derivs( &rxn, & reactions::diffeq );
        Jacobian     jacobn( &rxn, & reactions::jacobn );
        
        ////////////////////////////////////////////////////////////////////////
        //
        // Auxiliary timings
        //
        ////////////////////////////////////////////////////////////////////////
        eta          progress;
        if( tmax < dt ) tmax = dt;
        const size_t num_iter = Ceil(tmax/dt);
        
        const size_t every    =  (num_iter * dt_save) / tmax;
        
        ////////////////////////////////////////////////////////////////////////
        //
        //
        //                        SIMULATIONS
        //
        //
        ////////////////////////////////////////////////////////////////////////
        
        if( true )
        {
            ////////////////////////////////////////////////////////////////////
            //
            // Explicit constrained scheme
            //
            ////////////////////////////////////////////////////////////////////
            std::cerr << "-- Initializing Explicit Constrained Scheme" << std::endl;
            initialize(W); 
            save_fields( "chem-ini.dat", W);
            save_db( "chem-db.gz", 0.0, W,false,-1);
            num_newton = 0;
            const double newt_factor = 1.0/(volumes-1);
            
            
#if GFX == 1
            win->redraw(); Fl::check();
#endif
            std::cerr.flush();
            progress.reset();
            show_time(0,progress,0);
            
            {
                ios::ocstream infos( "chem-infos.dat", false );
                infos("#percent t_diff newt_ratio error\n");
            }
            
            double t_diff = 0;
            for( size_t iter=1; iter <= num_iter; ++iter )
            {
                const Real t     = iter * dt;
                const double td  = explicit_constrained_diffusion(W, dt);
                t_diff  += td;
                if( 0 == (iter%every) )
                {
                    const double frac_done = double(iter)/num_iter;
                    const double t_cycle   = t_diff/every;
                    progress(frac_done);
                    show_time(t,progress,1/t_cycle);
                    const double stamp = t;
                    save_db( "chem-db.gz", stamp, W,true,t_cycle);
                    {
                        ios::ocstream infos( "chem-infos.dat", true );
                        infos("%g %g %g %g\n", stamp, t_cycle, (num_newton*newt_factor)/iter, find_max_error(W));
                    }
                    t_diff = 0;
#if GFX == 1
                    win->redraw();
                    Fl::check();
#endif
                }
                
            }
            
            progress(1);
            const double t_end =  num_iter * dt;
            show_time( t_end, progress, 0 );
            fprintf( stderr, "\n"); fflush(stderr);
            save_fields( "chem-end.dat", W);
            save_db( "chem-db.gz", t_end, W,true,-1);
#if 0
            const double NewtRatio = (num_newton*newt_factor)/num_iter;
            {
                ios::ocstream infos( "chem-infos.dat", true );
                infos("%g %g %g %g\n",t_end, t_diff,NewtRatio,find_max_error(W));
            }
#endif
            
        }
        
        if( true )
        {
            ////////////////////////////////////////////////////////////////////
            //
            // Explicit standard scheme
            //
            ////////////////////////////////////////////////////////////////////
            double t_diff = 0;
            double t_reac = 0;
            
            std::cerr << "-- Initializing Explicit Standard Scheme" << std::endl;
            initialize(W);  
            save_fields( "norm-ini.dat", W);
            save_db( "norm-db.gz", 0.0, W,false,-1);
            
#if GFX == 1
            win->redraw();
            Fl::check();
#endif
            std::cerr.flush();
            progress.reset();
            show_time(0,progress,0);
            
            {
                ios::ocstream infos( "norm-infos.dat", false );
                infos("#percent t_diff t_reac t_all error\n");
            }
            
            for( size_t iter=1; iter <= num_iter; ++iter )
            {
                const Real t_old = (iter-1) * dt;
                const Real t     = iter * dt;
                const Real td    = explicit_standard_diffusion(W,dt);
                t_diff += td;
                double tr = 0;
                if( explicit_reaction )
                {
                    tr = reaction(W,t_old,dt,odeint,derivs,Y);                
                }
                else 
                {                    
                    tr = stiff_reaction(W,t_old,dt,stiff_odeint,derivs,jacobn,Y);
                } 
                
                t_reac += tr;
                const double t_cycle = t_reac + t_diff;
                
                if( 0 == (iter%every) )
                {
                    const double frac_done = double(iter)/num_iter;
                    progress(  frac_done );
                    show_time(t,progress,every/t_cycle);
                    ios::ocstream infos( "norm-infos.dat", true );
                    const double stamp = t;
                    infos("%g %g %g %g %g\n", stamp, t_diff/every, t_reac/every, t_cycle/every, find_max_error(W) );
                    
                    save_db( "norm-db.gz", stamp, W,true,t_cycle/every);
                    t_reac = 0;
                    t_diff = 0;
#if GFX == 1
                    win->redraw();
                    Fl::check();
#endif
                }
                
                
                
            }
            
            
            
            progress(1);
            const double t_end = num_iter * dt;
            show_time( t_end, progress,0 );
            fprintf( stderr, "\n"); fflush(stderr);
            save_fields( "norm-end.dat", W);
            save_db( "norm-db.gz", t_end, W,true,-1);
#if 0
            {
                ios::ocstream infos( "norm-infos.dat", true );
                infos("%g %g %g %g %g\n", t_end, tmx_diff, tmx_reac, tmx_diff+tmx_reac, find_max_error(W) );
            }
#endif
            
        }
        
        ////////////////////////////////////////////////////////////////////////
        //
        //
        //         POST PROCESSING  
        //
        //
        ////////////////////////////////////////////////////////////////////////
        
        compare_db();
        
        return 0;
    } 
    catch( const exception &e )
    {
        std::cerr << e.what() << std::endl;
        std::cerr << e.when() << std::endl;
    }
    catch (...) 
    {
        std::cerr << "Exception!" << std::endl;
    }
    return -1;
}